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Abstract 

In the substitutionally disordered narrow-gap 

semiconductor Pbi_2^^Gea;Te, a finite-temperature cubic-rhomb ohedral transi- 
tion appears above a critical concentration x ~ 0.005. As a first step towards 
a first-principles investigation of this transition in the disordered system, a 
(hypothetical) ordered cubic Pb3GeTe4 supercell is studied. First principles 
density-functional calculations of total energies and linear response functions 
are performed using the conjugate-gradients method with ab initio pseudopo- 
tentials and a plane-wave basis set. Unstable modes in Pb3GeTe4 are found, 
dominated by off-centering of the Ge ions coupled with displacements of their 
neighboring Te ions. A model Hamiltonian for this system is constructed 
using the lattice Wannier function formalism. The parameters for this Hamil- 
tonian are determined from first principles. The equilibrium thermodynam- 
ics of the model system is studied via Metropolis Monte Carlo simulations. 
The calculated transition temperature, Tc, is approximately 620K for the cu- 
bic Pb3GeTe4 model, compared to the experimental value of Tc ~ 350/'^ for 
disordered Pbo.75Geo.25Te. Generalization of this analysis to the disordered 
Pbi_a;Gea;Te system is discussed. 
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I. INTRODUCTION 



Numerous substitutionally disordered ferroelectrics exhibit structural phase transitions 
in which the nature of the phase transition changes as a function of composition. Examples 
include Bai_^.Sr^.TiO30, PbZri_^Ti^.03il, and Pbi_^Ge^Tel Pbi-^Ge^^Te is an ideal system 
for the development of first-principles methods to investigate the effect of substitution on 
structural phase transitions for several reasons. First, the unit cells of the endpoint com- 
pounds are small. Pure PbTe and GeTe have only two atoms and ten valence electrons per 
unit cell, as opposed to five atoms and 24 or more valence electrons per unit cell in the per- 
ovskite oxide ferroelectrics. Secondly, the phase diagram in Pbi-^^GcxTe is rather simple.0 
Pure PbTe has a rocksalt structure that is stable to zero temperature. Pure GeTe also 
has the rocksalt structure at high temperature but undergoes a structural phase transition, 
indicated in Figure [^, to a rhombohedral phase at a critical temperature ~ 670K.i As 
the Ge concentration x decreases from 1 to 0, for a cubic-rhombohedral phase transition 
decreases smoothly to zero at x ~ O.OOS.Si Finally, there are existing phenomenological and 
empirical models for this transition.BHi^ These models focus on the role of Ge off-centering0'i 
in the phase transition. Off-centering is observed in a number of compounds where there is 
a mismatch between the ionic radii of two species that statistically occupy the same kind of 
site.0 There is direct experimental evidence for off-centering in Pbi_2;Ge2;Te. The extended 
X-ray-absorption fine-structure measurements of Islam and BunkerS show two peaks in the 
distribution function for Ge-Te distances both above and below Tc, as opposed to the single 
peak that would be seen if the Ge atoms were located at the centers of the octahedra formed 
by their first-neighbor Te atoms. The phenomenological and empirical models previously 
considered provide a reference for comparison with the model derived from first principles 
here. 

Model Hamiltonians based on ab-initio calculations have been successfully constructed for 
a number of stoichiometric ferroelectrics and related materials, including GeTfl PbTioB, 
PbZrOaB, BaTiOsliE, SrTiOs^ and KNbOa^. In constructing these models, the high- 
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symmetry average structure of the high-temperature phase is chosen as a reference struc- 
ture. Normal mode dispersion relations and eigenf unctions are calculated to determine the 
unstable modes. The normal mode branches containing the unstable modes form a basis for 
the ionic displacement subspace determining, by projection, an effective Hamiltonian. Via 
a linear transformation, a localized lattice Wannier function^ (LWF) basis is found. The 
energy is expanded in powers of the lattice Wannier function coordinates and strain and the 
coefficients in the expansion are determined from a set of ab-initio calculations. 

In a substitutionally disordered system, it is not enough to include only the lattice degrees 
of freedom. Configurational entropy also plays an important role; the partition function, 
from which all thermodynamic properties can be obtained, includes a sum over all possible 
configurations. Different fixed configurations will undergo structural phase transitions at 
different temperatures. Furthermore, the nature of the configurations whose properties 
dominate in the thermodynamic limit will itself depend on temperature. However, it is not 
necessary to investigate all possible configurations in order to model a phase transition in a 
disordered system. In fact, if a model can be formulated in terms of intersite interactions that 
decay rapidly with distance to some asymptotic form, only a small number of configurations 
need to be explicitly calculated from first principles. 

In this work, we investigate a single Pbi_2,GexTe configuration: the ordered 8 atom 
Pb3GeTe4 cubic cell {x = 0.25), shown in Figure This configuration was chosen for 
several reasons: (1) it is toward the Pb-rich end of the phase diagram, which is appropriate 
for studying how adding a small concentration of Ge to PbTe "turns on" a phase transition 
at zero temperature, (2) it contains the minimum number of atoms per cell (eight) for 
any x < 0.25 and (3) among the eight atom cells with x = 0.25, it has the highest point 
symmetry. The small size and high symmetry of the unit cell for cubic Pb3GeTe4 makes 
ab initio calculations relatively computationally inexpensive. Our calculation of for this 
system provides a benchmark for comparison with Tc for other configurations and for the 
ensemble average at the same composition. 

This paper is organized as follows. Section II describes the general principles governing 
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the construction of model Hamiltonians and the lattice Wannier function method for deter- 
mining the effective Hamiltonian subspace. In Section III, these methods are applied to the 
specific case of cubic Pb3GeTe4 and a model Hamiltonian is obtained. The model Hamil- 
tonian is used in a classical Monte Carlo simulation in Section IV to obtain the transition 
temperature and order of the cubic- rhombohedral phase transition in this system. The re- 
sults are further discussed in Section V, including their implications for the substitutionally 
disordered Pbi-^GeajTe system. Conclusions are given in Section VI. 



II. CONSTRUCTION OF MODEL HAMILTONIANS 

At the level of the Born-Oppenheimer approximation for electronic energy, the classical 
partition function Z for a disordered system on a fixed lattice is 

Z=Y.J (^{^^j} J / de exp[{-(3E{{aj}, {u,}, {ii,}, e)], (2.1) 

where {aj} represents the chemical configuration, i.e. which type of ion occupies each site 
j, e is the homogeneous strain tensor, {uj} is the set of ionic displacements defined with 
respect to some reference structure for the given chemical configuration and strain, and 
{iij} is the set of ionic velocities. If, in addition, the ionic motion is treated classically, the 
integral over {uj} leads to a trivial u-independent factor, and the partition function can be 
rewritten as 

I de exp[(-/5E({a,}, {u,}, e)], (2.2) 



The partition function |2.2| can be rewritten Z oc J2{aj} ZdcTj}), where Z{{aj}) is the 
partition function for the ensemble corresponding to the single chemical configuration {crj}- 
In what follows, we consider only one chemical configuration, so the configuration label aj 
will be dropped. The partition function now simplifies to 

Z J d{uj} J de exp[{-PE{{uj},e)]. (2.3) 
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In general, E{{uj},e) can be expanded as a Taylor series in powers of uj and e. If the 
reference structure is at an energy extremum, then the linear terms in the expansion vanish: 

^({u,}, e) = + E'^'\{u,}, e) + E'^'\{u,}, e) + . . . . (2.4) 

The harmonic term, E'^'^\{uj}, e) is the sum of a lattice term 

Fijap Uia Ujf3, (2.5) 

ija0 

where F is the force constant matrix, a strain term 

C„^7<5 Gap e^s, (2.6) 
where c is the elastic constant tensor, and a strain-coupling term 

diaPy Uia ^13^, (2.7) 

io/37 

where g is the strain-coupling tensor. By the usual change of variable to normal mode 
variable, {u^} = J^u^^u^uy the harmonic lattice term is reduced to a single sum: 

FijafB Uia Uji3 = ^ ^m,,t^^a^, (2.8) 

ijafi ^ 

where ay is the amplitude, ey the (normalized) ionic displacement pattern, ujy the frequency, 
and rriy the mode effective mass for normal mode v. 

In what follows, it is assumed that the fixed configuration is periodic. Then the normal 
mode label v can be replaced by kia, where k is the wavevector, i the symmetry-invariant 
subspace and a the branch of the normal mode dispersion curves on which the mode lies. 
For example, in a structure with two atoms per unit cell, i = 1 represents the acoustic 
modes, i = 2 the optical modes, and in each subspace a = 1 labels the longitudinal and 
a = 2,3 labels the two transverse branches. 

If the extremum or reference structure is a saddle point in the energy surface rather than 
a minimum, some of the normal modes will be unstable (c^kia < 0). As described in detail by 
Rabe and WaghmareHil, the model Hamiltonian approach to structural phase transitions 



associated with soft phonons aims to reduce the ionic degrees of freedom to those of the 
"effective Hamiltonian" subspace that contains the important anharmonic terms in Eq. |2.4 
All normal mode subspaces i that contain unstable modes must be included in the effective 
Hamiltonian subspace. An expansion of the total energy in these degrees in freedom requires 
higher-order terms in order to stabilize the structure at a finite distortion. As the simplest 
approximation, terms are kept to only harmonic order in the subspace complementary to the 
effective Hamiltonian subspace and higher order mixing of the effective Hamiltonian subspace 
with the complementary subspace is neglected. This approximation is valid to the extent 
that the magnitude of the neglected higher order terms is small for all ionic displacement 
patterns that contribute significantly to Z, i.e. those where E{{uj},e) is small. Integration 
over the complementary subspace gives a structure-independent contribution to Z which 
can be neglected. The remaining terms in E give an effective Hamiltonian 'Heff- Integration 
of Ti-eff over the effective Hamiltonian subspace is sufficient to reproduce the dependence of 
the structure on temperature. 

A change in variables allows the basis of the effective Hamiltonian subspace to be con- 
verted into one where the displacements are localized, i.e a. lattice Wannier function (LWF) 
basis wrEI. This procedure is analogous to that which des Cloizeauxill gave for determin- 
ing electronic Wannier functions for a multidimensional subspace; except that the periodic 
functions are now ionic displacement patterns.0 For simplicity, we assume that the effective 
Hamiltonian subspace consists of only a single normal mode subspace « = 1, drop the sub- 
script i, and let label the normal modes in this subspace. The LWF is given in terms of 
the normal modes by 

= /^^rfk(5:a^(k)6k.e^('^--'^-^)). (2.9) 

The subscript (3 labels the components of the LWF. The number of components of w/j is 
equal to the number of branches in the effective Hamiltonian normal mode subspace. The 
matrix C (k) is included so that the polarizations of the different normal mode branches are 
correctly related to the LWF components across the Brillouin zone. The phases 0ka are 

6 



very important. They are chosen to make the ionic displacement pattern corresponding to 
w both as local as possible and with as high a symmetry as possible. By construction, the 
LWF has translational symmetry: Tr,/(w/3r) = tf/jR+R/, where R' is any lattice vector. 
The normal modes are given in terms of the LWF basis via the inverse transform 

[e^^'-^Wa = Y.^C-Xpw^^e''''^. (2.10) 

R/3 

Any displacement pattern Uj in the effective Hamiltonian subspace can be written as 

Ui = H'^kaeka- (2.11) 

Using Eq. |2.1CI| , and changing the order of the summation, we obtain 



u,- = E ^/3R[E(^"')«/5«kae^(^-^-'^'''^)] = E ^0^i0n. (2.12) 

R/3 ka R/3 



The new vector parameter ^r gives the displacement field via Eq. p. 12 and provides the 
basis for a spin-like model of the system. 

III. MODEL HAMILTONIAN CONSTRUCTION FOR PB3GETE4 

In this section, we describe the construction of the effective Hamiltonian subspace and 
the explicit expression for the model Hamiltonian for Pb3GeTe4 from first-principles density- 
functional calculations within the local-density approximation (LDA). Bachelet, Hammann 
and Schliiter pseudopotentialsii were used for each species in the two-projector Kleinman- 
Bylander form.il@ For Pb and Te, the d potential was taken as local; for Ge, a d local 
potential was found to give spurious low-energy "ghost" stateJ^lii and so a p local po- 
tential was chosen. Total energy calculations were done using the CASTEP 2.1 program of 
Payne et aM CASTEP 2.1 performs conjugate gradients minimizationil of the LDA density 
functional total energy in the Kohn-Sham formalism.^ The exchange-correlation energy was 
Perdew and Zunger's parameterizatiorS of the Ceperley-Alder valued^ for the uniform elec- 
tron gas. The calculations were performed for a 8-atom unit cell with a 300 eV cutoff for the 
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plane wave basis set (about 3000 plane waves per k-point), a 4 by 4 by 4 Monkliorst-Pac 
k-point set, and a 36 by 36 by 36 real-space grid for the charge density. Force constant matri- 
ces at the high-symmetry q 7^ points in the Brillouin zone and Born effective charges were 
computed using density-functional perturbation theoryilS in the variational formulation.ii 
At q = 0, the force constant matrix was computed using Hellmann-Feynman forces, as de- 
scribed in more detail below. Pulay correctionsil were added to the total energy results as 
described in the Appendix. 

We began with Pb3GeTe4 in a high-symmetry reference structure. The structure, shown 
in Figure |], is produced by replacing a cubic superlattice of Pb ions in the PbTe rocksalt 
structure with Ge ions. There are eight atoms per unit cell. The structure has full cubic 
symmetry (space group Pm3m). The Ge occupies the cell corner (Wyckoff position 1(a)) 
and the three Pb's occupy the face centers (3(c)). The four Te's occupy two crystallograph- 
ically distinct positions: the edge centers (3(d)) and the cube center (1(b)). Although this 
structure is not the minimum-energy structure, the forces on all ions are zero by symmetry, 
and the problem of relaxation can thus be neglected (note this is not true for the general 
case of partial Pb Ge substitution). Furthermore, the Pb3GeTe4 structure with the lat- 
tice constant a which minimizes the total energy is at an extremum of the E{{ui}, e) energy 
surface, so that it can be used as a reference structure in the analysis of Section II. We found 
a = 6.275 A for cubic Pb3GeTe4. Although there is no experimental cubic Pb3GeTe4 sample 
with which this result can be prepared, we can estimate the "experimental" lattice parameter 
in two different ways as follows. First, we extrapolate the experimental lattice parameters 
for PbTe at room temperature0 and for GeTe at the cubic-rhombohedral phase transition 
temperature^ to zero temperature via the thermal expansion coefficient of PbTe at room 
temperature^ and then average the results assuming Vegard's law to obtain a = 6.307(2) A. 
Second, we estimate the experimental lattice parameter for disordered Pbo.75Geo.25Te at 
room temperature and the thermal expansion coefficient for nearby compositions via the 
figures in Reference ^ and then extrapolate to zero temperature to obtain a = 6.311(1) A. 
The two estimates are in good agreement. The LDA lattice parameter is about 0.5% lower 
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than the "experimental" value. This magnitude of lattice constant underestimate is typical 
of LDA calculations.!! 

Next, we identified the lattice instabilities of the reference structure. The force constants 
at q = were found by displacing each ion in turn and then calculating the Hellmann- 
Feynman forces. Anharmonic effects were found to be significant in Pb3GeTe4 and thus 
it was necessary to calculate the forces for small (0.001 to 0.01 A) displacements and ex- 
trapolate the force constant values to zero displacement in order to obtain high precision. 
For comparison, we also calculated the force constant matrix at q = using the density- 
functional perturbation theory method. The latter results, however, violated the acoustic 
sum rule by as much as 0.1 eV/A^, whereas the maximum violation of the acoustic sum 
rule for the extrapolated Hellmann-Feynman forces was only lO"'' eV/A^. The discrepancy 
between the force constant matrix elements found by the two methods was largest for the 
diagonal terms. In order to avoid significant acoustic sum rule violation corrections, at q = 
we used the extrapolated Hellmann-Feynman forces. 

Diagonalization of the resulting q = dynamical matrix yielded a single unstable mode 
at F (see Table |), with symmetry F15. The fact that it is strongly dominated by Ge motion 
led us to choose a single LWF centered on Ge with F15 (vector) symmetry. To construct 
the effective Hamiltonian subspace, it is therefore necessary only to investigate those q 7^ 
normal modes compatible with this choice of LWF symmetry. We performed q 7^ linear- 
response calculations to calculate the normal mode frequencies and eigenvectors for all such 
modes at the BZ points F, R, X and M (Table I). For each label§ X[,Xi„ M[, M'^ and 
i?i5, the lowest energy mode also has a large component of Ge motion (see Tables |I|and 0), 
and is therefore easily identified as a mode belonging to the effective Hamiltonian subspace. 

Substitution of the first-principles results for the ionic displacements of the selected nor- 
mal modes into Eq. |2.10| leads to a set of linear equations, which can be solved to obtain the 



LWF. To obtain the exact LWF, which in general involves nonzero ionic displacements out to 
infinite distance, it would be necessary to include an infinite number of independent normal 
modes in the analysis. We expect, however, that the ionic displacements corresponding to 
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the lattice Wannier function decrease rapidly with distance, in analogy to the behavior of 
electronic Wannier functionsS and of the lattice Wannier function for PbTiOa.El. A good 
approximation to such a LWF can be obtained by setting all ionic displacements equal to 
zero outside of some finite region. The remaining finite set of ionic displacements can be 
determined from a finite set of normal modes via Eq. ^.9| . For cubic Pb3GeTe4, we used all 
the modes of Table ||. All modes were normalized so that the sum of the squared ionic dis- 
placements was 1 A^ per primitive cell. For these high symmetry points, the normal modes 
can be chosen so that Ge motion is always strictly along a Cartesian direction and thus the 
matrix C in Eq. p.9| is the unit matrix. The phases 0ka are already incorporated into Table 
m and are chosen so that Ge motion is always in the same direction. The included modes 
yielded a set of 21 independent linear equations and therefore allowed the 21 independent 



components of ionic displacement shown in Table |T| to be determined. The displacements 
corresponding to the approximate LWF involve all ions out to {V5/2)a from the central Ge 
ion and some ions as far as 1.5 a away. 

As can be seen from Table PTT| and Figure |^, the LWF satisfies the locality assumption very 
well, with first-near neighbor Te displacements only 0.27 times the central Ge displacement, 
second near- neighbor Pb displacements up to 0.14 times the Ge displacement and all further- 
neighbor displacements less than 0.05 times the central Ge displacement. Note that although 
the LWF component shown transforms as the z component of a vector, individual ionic 
displacements, such as the Pb displacements shown in Figure can be in directions other 
than z. Henceforth, we use units in which a dimensionless local distortion amplitude of 1 in 



the z direction corresponds to the displacements shown in Table |lTl . 

To summarize, the structure-dependent part of the partition function for Pb3GeTe4 
should depend only on a subspace of the full ionic displacement space, which can be used to 
generate an effective Hamiltonian. So far, we have found a good approximation to wr, the 
lattice Wannier function basis for the effective Hamiltonian subspace, which includes the 
unstable normal modes. We have a spin-like representation in which there is a one-to-one 
correspondence between ionic displacement patterns in the effective Hamiltonian subspace 
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and the magnitudes and directions of a set of vectors located on each Ge ion. We will 
now describe the model Hamiltonian for Pb3GeTe4, obtained by (1) expanding the effective 
Hamiltonian energy per unit cell 7ie//({^i}, e)/A^ in powers of C,i and strain e, (2) truncat- 
ing this expansion, and (3) determining the coefficient for each term from first-principles 
calculations. 

Previous ab initio studies on perovskite ferroelectrics such as BaTiOslii and PbTiOslii 
have led to models of the same nature as that presented here, namely a set of interacting 
vector spins sitting on Wyckoff positions of full cubic symmetry for a cubic lattice with 
underlying space group Pm3m. For Pb3GeTe4, we used the same truncations in the ex- 
pansion of 'Heff as was used for the PbTiOa modehlli We verified that the terms retained 
were necessary and sufficient to accurately fit the corresponding ab initio results. Various 
contributions to the model Hamiltonian will now be considered in turn. 

The total energy contains a structure-independent constant term 

Uo/N (3.1) 

— * 

giving the energy of the high symmetry Pb3GeTe4 reference structure (^j = e = 0). In 
the case of disordered systems, differences in reference structure energy are important for 
comparing different chemical configurations. The value oIUq/N for cubic Pb3GeTe4 is given 
in the Appendix. Since only one cation configuration is considered in this work, the value 
of Uq/N does not affect subsequent results. 

The lowest order terms which depend on the local distortions {^j} are harmonic, with a 
local contribution Uih = Z]j^|6P and an intersite pair contribution Uph = CLijais^ia^jis- 
At long range, we expect the intersite pair interaction to approach the dipole-dipole form, 
while at shorter range, there should be significant corrections due to higher order multipole 
effects, induced charge redistributions and direct overlap of the ionic displacement patterns. 
We thus split Uph into the sum of a long-range dipole-dipole term, Udd, and a short range 
"correction" term, Ugr- 

The energy of a system of dipoles of moment /I in a medium of electronic dielectric 
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constant 600 is given by:0 

UMfi.})/N 4 EE - (3.2) 

The sum over d is a sum over all interdipole separations, /Tj is the dipole moment at site z, 
and /ij_,_^is the dipole moment of the dipole at distance d from the dipole at site i. 

To make use of Eq. p.2| , we need to know the dipole moments /ij corresponding to local 
distortions ^j. Since the dipole moment at ,^ = is zero via centrosymmetry, the dipole 
moment for nonzero ^ is simply the sum of the polarizations Pj induced by individual ionic 
displacements Uj specified by the given ^i. 

fi = J2Pji^M))- (3.3) 

j 

For small displacements, the Pj are determined by the Born effective charge tensor for the 
ion j, defined as the differential change in total polarization due to displacement of that ion: 

Symmetry considerations reduce the number of independent terms in the Born effective 
charge tensors for Pb3GeTe4. Though our Pb3GeTe4 structure has cubic symmetry, not all 
of the ions sit at centers of full cubic symmetry. The Ge and Te(l) ions sit at centers of cubic 
symmetry, while the Pb and Te(2) ions fill the two Wyckoff positions of multiplicity three, 
each of which is a center of tetragonal symmetry. The Born effective charge tensor is still 
diagonal under tetragonal symmetry, but the axial and perpendicular values are different. 
Thus, there are 6 independent Born effective charge tensor components in the Pb3GeTe4 cell, 
one each for Ge and Te(l) and two each for Pb and Te(2). The Born effective charge tensors, 
obtained from linear response calculations for cubic Pb3GeTe4 at the equilibrium LDA lattice 
constant, are shown in Table |V|. The sum of the Born effective charges + 
as determined via linear response calculations was —0.36. This violation of charge neutrality 
is a consequence of finite k-point sampling.ii We corrected for it@ by adding the same 
constant 0.045 to all diagonal Born effective charge tensor components. The diagonal cation 
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(Te) effective charge components are in the range +6 to +8 (—6 to —8). These effective 
charges are anomalously large compared to the formal valence charge states, but consistent 
in magnitude with those found experimentally^ and theoretically^ for binary group IV 
tellurides. 

Since the Uj are proportional to (Eq. ^.12| ), the dipole moment fl can be written simply 

as 

iI, = J2Z^n^i^^=Tii,), (3.5) 

j 

where symmetry leads to a scalar mode effective charge parameter Z* , which for Pb3GeTe4 
has the value 12.10 eA. Eq. p.2| can then be written in terms of ^ as follows: 

u^am/N 4 E E (£2 (3.6, 



e 



The other material specific parameter needed in Eq. p.2| is eoo- By symmetry, the 
dielectric tensor is isotropic in cubic Pb3GeTe4. Via linear response methods, we calculated 
= 46.7. The experimental value for PbTe is 33-34 at room temperatureSS and 40 
at zero temperature.0 For GeTe, the measured valued at room temperature is 36. For 
Pbi_^.Ge^Te, eoo measurements are available for x = 0.07, for which the values obtained are 
about 39 at room temperature, 41 for the cubic phase just above the phase transition and 
40 at zero temperature.0 Our value for eoo is somewhat higher than the value obtained from 
interpolation of the experimental values, consistent with the general tendency of LDA to 
overestimate dielectric constants.^ 

Next, consider Usr, the short-range correction term in the intersite pair interaction. The 
force constant matrix elements determined in a series of q 7^ linear response calculations, 
along with the corresponding ionic displacement fields, allow the total harmonic energy 
Uih + Uph to be determined. These are given in Table |V| for the various high-symmetry q 
points used in determining the LWF. For each normal mode used to construct the LWF, 
we calculated the dipole-dipole interaction energy Udd via Eq. We then subtracted the 
dipole part from the total harmonic energy to obtain Usr{{^i})- 
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As in the case of the LWF itself, the information obtained from a finite set of independent 
q points hmits the number of interaction coefficients that can be determined. We include 
here all symmetry-adapted couplings to third neighbors in the cubic lattice (or sixth neighbor 
cations). There are two independent coefficients at first neighbor 

uAm/N = ^ E E [^^te ■ d){i+i ■ d) + &Tte • - ■ ^kC,- ■ m, (3.7) 

* d=nnl 

three at second neighbor 

+ ^ E E \dL{i ■ ■ d) + dTi{^ ■ d,){i^^ ■ d,) 

* d=nn2 

+ c?T2(6-4)(5+r4)], (3.8) 
and two at third neighbor 

+ ^ E E [/^te ■ mu ■ d) + /t(6 ■ - • d){iu ■ d))]. (3.9) 

* d=nn'i 

In these expressions, c? is a unit vector that is summed over all directions from a site i to 
the neighbors in a given shell. For second neighbors, a further distinction is made between 
the two direction orthogonal to d: d2 is the Cartesian vector orthogonal to d, and di is 
orthogonal to both d and ^2- ii^d is the value of ^ for the site at distance cT from site i. 
By substituting the results of Table and the appropriate values of into Eqs. |3l7| - |3l9| . 



we obtained the local harmonic coefficient A and the intersite coupling constants 6l, 6t, 



di + dxii dj'2 and /l + "^fr shown in Table VI. There was not enough information at the 



high symmetry q points used to separate d^ from dxi or from fx- As long as di, dxi, 
fi and fx were all nonnegative, however, we found that the values chosen had little effect 
on the thermodynamic properties of our model, so we set = fr and d^ = dn- Note that 
the local correction to the dipole-dipole energy decreases as the intersite distance increases, 
becoming smaller than the dipole-dipole interaction at third neighbors. While there is a 
tendency for the intersite interaction to approach the dipole-dipole form at large distances, 
there is still significant deviation from this form at third neighbor distance. 
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Next, consider the anharmonic terms to be included in the model Hamiltonian. Such 
terms are crucial for determining the thermodynamics of a system undergoing a mode- 
softening structural phase transition. Making an important simplifying approximation which 
has been previously applied in the study of GeTeS and perovskite oxides^^S, we ignore 
multisite and higher-order pair interactions in our model Hamiltonian and include anhar- 
monic terms only in the local distortion energy Uia- The Ge sites are centers of full cubic 
symmetry; thus the onsite energy expansion includes only cubic invariant terms. We trun- 
cate the expansion at eighth order, including all terms up to quartic order and the isotropic 
sixth and eighth order terms: 

Uia/N = ^ Y.im' + + ity + ii) + + EU'). (3.10) 

The coefficients B, C, D and E were obtained from a series of frozen phonon calculations 
at r, after subtracting the V point harmonic energy, — 1.6524 eV (Table 0). Various 
amplitudes of distortion |,^| were applied in both the z and {x + y + -2)/a/3 directions and 
a least squares fit to Eq. |3.10| was performed using the ah initio total energy results. The 



resulting coefficients are given in Table 0. A comparison of the frozen phonon and fitted 
results is shown in Figure |5(a). A contour plot of the equal energy surfaces of the fit is 
shown in Figure §(b). The energy minima are located along the [111] directions. 

While the spin-like model obtained so far includes the necessary ingredients for a phase 
transition, the ground state symmetry of the crystal, degree of distortion and the nature 
of the phase transition all depend on magnitudes of both the strain (Eq. |2.6| ) and strain 
coupling terms (Eq. ^.T]).!^ We include the effects of strain and strain couphng to lowest 
order. The symmetry adapted forms are: 

Ustrain/N = ^Cu + l^u ^aa^PP + 7^44 Y ^1/3 (S-H) 



and 



U strain coupling / N 7y (^I] ^aa) l^i I lyr 'YX^cea ^jq) 



+ I7 E (3-12) 
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The variations in cell energy of uniformly strained Pb3GeTe4 structures, with Pulay 
corrections, are shown in Figure ^ for the following types of strain: uniform dilation [e^x = 
€-yy = Gzz), tctragoual strain {cxx = ^yy = ^ ^zz/'^) and rhombohedral shear {cxy = Cyx = Cxz 
= Gzx = = Czy). Least squares fits to these results determine the elastic constants Cu, 
Ci2 and C44 shown in Table |VI[ 

Finally, strain couplings were calculated by subtracting the energy at ^ = from various 



uniform fixed amplitude (|^| = 0.25) distortions under various strains, using Eq. |3.12 
to determine the form of the corresponding energy differences and finding the unknown 
coefficients via the linear parts of least-squares quadratic fits to each curve. The results are 



shown in Fig. In Table |V|, we present the values obtained for the three independent 
strain coupling parameters g^, gi and (72- As a check, further calculations were performed 
in which the amplitude of the LWF distortion was also varied. The results were consistent 
with the strain coupling parameters found. 

The effect of strain on the energy of a F local polar distortion is shown in Figure |^, which 
is similar to Figure ^(b), except that the distortion energy is now minimized with respect to 
uniform strain. There are still 8 local minima along the [111] directions. Their energies are 
17 meV/cell lower, and corresponding values of 0.035 larger than if the terms involving 
strain are neglected. 

IV. FINITE TEMPERATURE SIMULATIONS 

The model Hamiltonian constructed in the previous section is completely specified by 

U ({Ci}, e)/A^ = {Uo + Ulh + Udd + Usr + Ula + Ustrain + Ustrain coupling) / N , (4.1) 



where the individual contributions are given by Eqs. |3.1| and |3.6[ - p.l2| and the parameters 



appear in Table [V|. This model Hamiltonian applies to all ionic displacements believed 
to play a role in the structural phase transition and allows one to calculate Tc and other 
properties in the vicinity of the transition for the ordered Pb3GeTe4 system. We analyzed the 
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finite-temperature behavior using a classical "single flip" Metropolis Monte Carlo simulation. 
Three types of Monte Carlo steps were applied: (1) A "spin" was chosen at random. A 
random vector located within a cube centered on the origin was added to the chosen spin 
vector. Acceptance ratios of about 1/e were obtained for a cube radius of 2VkT, where kT 
is in eV. (2) A spin was chosen at random. The x, y and z components were independently 
multiplied by ±1 with probability 0.5. This led to a great increase in simulation speed as the 
spins were able to overcome barriers between local energy minima. An acceptance ratio of 
about 0.3 was found for this step near Tc. (3) One strain component was chosen at random. 
A random number between — e and +e was added to this strain component. Acceptance 



ratios of about 1/e were obtained for cubic cells for e = O.lSyfcT/L, where kT was in eV 
and L was the size of the cube edge in units of a. 

We simulated an L x L x L supercell of cubic Pb3GeTe4, with L = 10, representing 
8000 atoms. We determined mean total energies and lattice parameters as a function of 
temperature. At each temperature, the parameters were averaged over A^^ Monte Carlo 
attempts for Nc = 100 cycles. The longest autocorrelation time for the q = (27r/L)z Fourier 
component of the spin fleld was determined to estimate the decorrelation time r at each 
temperature. For a proper ensemble average, NgNc >> r. We generally selected A^^ so that 
NgNc ~ 20r. From the root mean square deviation of the individual cycle averages, arms, 
we can obtain an estimate of the rms error for the ensemble averaged quantitiesS, 



We ran the simulations for a range of temperatures between OK and lOOOK. We began at 
lOOOK, cooled to OK and heated back up to test for hysteresis effects. At each temperature, 
the system was allowed to equilibrate before any sampling began. 

At each temperature, the symmetry of the average cell was either cubic or rhombohedral 
to within statistical errors. Therefore the lattice parameter and the rhombohedral angle as a 
function of temperature are sufficient to give the unit cell as a function of temperature. The 
results are shown in Figure || and are consistent with a second order cubic-rhombohedral 





(4.2) 
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phase transition at Tc = 620K. The error bars for certain temperatures very near Tc may 
be consistent with shght hysteresis, indicating a very weak first order transition. However, 
since this is observed in the temperature range where the effects of the periodic boundary 
conditions are significant, no firm conclusion can be reached. In any case, the hysteresis 
effect, if it exists, spans less than 10 degrees. Experiments on disordered Pbi„a;Gea;Te for 
small X show a continuous transitionH and give Tc ~ 350/^ for x = 0.25. This apparent 
discrepancy in Tc will be discussed in detail in the next section. 

The plot of lattice parameter versus temperature shows a discontinuous slope at Tc 
and a negative coefficient of thermal expansion below Tc. The model cannot be expected to 
obtain the right value for thermal expansion coefficient because it neglects anharmonic terms 
involving the normal modes complementary to the effective Hamiltonian subspace. The 
inclusion of additional modes should not alter the predicted discontinuity in the expansion 
coefficient, however, and a discontinuity of the right magnitude has indeed been observed 
experimentally in Pbi_3;.Gea,Te.@ 

V. DISCUSSION 

In this section, we discuss further various features of our model Hamiltonian, compare 
the model with other systems and models, comment on the apparent discrepancy between 
our Tc and the experimental value, and discuss what the present calculations have taught 
us about generalizing to the case of disordered systems. 

The harmonic part of our model is completely given by the normal mode dispersion 
curves. In Figure ^, these dispersion curves are shown for high symmetry paths in the 
Brillouin zone. Three features are particularly noteworthy: (1) Roughly speaking, there is 
a relatively fiat branch at u; ~ 425 cm~^ and a relatively fiat branch at u; ~ 600^ cm~^. 
This refiects the dominance of the longitudinal first neighbor dipole-dipole interaction and 
correction term 6^. (2) Although weaker than the first neighbor longitudinal interaction, 
the other interactions are crucial for determining the ground state structure. For example, 
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if the correction to the first neighbor transverse interaction, bx, were zero, the most unstable 
modes would have symmetry My and the ground state would be antiferroelectric rather than 
ferroelectric. (3) The long-range dipole-dipole interaction leads to LO-TO splitting at F. The 
frequency of the TO mode at T is still negative, as opposed to the case of various perovskite 
ferroelectrics@, where the LO mode corresponding to the soft TO mode generally has very 
high frequency. This is a consequence of the strength of the Ge off-centering instability in 
Pbi_a.Gea;Te. 

The form of our model Hamiltonian is that of interacting vectors on a cubic lattice. The 
same form of model Hamiltonian applies to perovskite ferroelectrics; thus the difference in 
phase transformation behavior of the different models can be related to the differences in 
the parameters and/or the number of vectors per unit cell. A comparative study of all the 
models to date is outside the scope of this paper; we will just compare the effects of strain 
coupling in our Pb3GeTe4 model and the PbTiOs model of Waghmare and Rabe. Ill The 
effect of strain coupling in Pb3GeTe4 is much weaker than in PbTiOs. For each case, we 
compared the energy in the ground state of the corresponding model with the symmetric 
reference state energy. In Pb3GeTe4, the distortion energy was -121 meV/cell, while strain 
and strain coupling contributed only -20 meV/cell. In PbTiOs on the other hand, the strain 
and strain coupling energy was -233 meV/cell at the ground state. Neglecting these terms, 
the distorted ground state was 120 meV/cell higher in energy than the reference state. The 
difference in the strength of strain coupling also makes the difference between a continuous 
phase transition in our Pb3GeTe4 model and a first-order phase transition in the PbTiOs 
model. 

In the Ge off-centering model of Pbi-^Ge^. each Ge atom is displaced off-center in 
one of the 8 cubic [111] directions. The cubic-rhombohedral phase transition is an order- 
disorder transition involving the directions of the displacements. Above Tc, the directions of 
the Ge atom displacements is disordered. There is no net polarization and the symmetry is 
cubic. Below Tc, the directions of the Ge atom displacements are ordered and lie predomi- 
nantly along one axis, lowering the symmetry from cubic to rhombohedral and making the 
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system polar. 

Our model, derived from first principles, essentially confirms the Ge off-centering model. 
It contains a local polar distortion (the LWF) dominated by Ge motion, whose energy 
minima are in directions where the Ge are displaced along the [111] directions (Figure ^ 
in agreement with the off-centering model. To see if the phase transition in our model is 
order-disorder, we calculated the distribution of local distortion vectors at Tc - 50K and Tc 
+ 50K. The results are shown in Figure |10[ Above Tc, there are distinct peaks in all eight 
[111] directions, while below Tc, there is only a single peak in one of the [111] directions, 
confirming the order-disorder nature of the transformation. We note here that EXAFS 
measurements of Islam and Bunkei0 leave an ambiguity in the direction of displacements 
of the Ge atoms at low temperature. For a ground state with rhombohedral strain along 



the (111) direction, say, both (111) and (111) local polar distortions lead to the same Ge-Te 
near neighbor pair distribution functions. In our model, the interactions between the local 
polar distortions have been determined and they are found to lead unequivocally to strictly 
ferroelectric ordering of these distortions. 

EXAFS measurements^ show two peaks in the Ge-Te distribution function in 
Pbi-j^Gca-Te. Only a displacement of Ge along one of the [111] directions breaks the symme- 
try in a way that leads to two different Ge-Te near neighbor distances. Although the local 
polar distortion of our model involves motion of atoms other than the central Ge, the same 
result holds that only a local polar distortion along a [111] direction leads to exactly two 
different Ge-Te distances. However, in our model, both the fact that the ionic displacement 
patterns corresponding to local polar distortion on different Ge atoms can overlap and the 
fact that the local polar distortions do not lie strictly along the [111] directions (see Figure 
|T0|) makes it less certain that there is a two-peak Ge-Te distribution function. By transform- 
ing the {C,i} coordinates of our model back to ionic displacements via equation |2.12 , 



we can 



calculate the first-neighbor Ge-Te distribution functions of our model. We show our results 
in Figure ITU|. There are two peaks, even above Tc, although there is significant overlap. 



These distribution functions do not incorporate the thermal noise of the neglected normal 
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modes, which would tend to further broaden the peaks; nonetheless, they are consistent with 
the experimental results showing two peaks. 

Since the semi-empirical models of Katayama and Muraseli (KM) and of Yaraneri et a/J 
are based upon Ge ions tunneling between (111) and (111) positions, in apparent contradic- 
tion with our results showing one preferred position below Tc, we have investigated this issue 
more carefully. Extrapolating our model to lower Ge concentrations by setting the dipole- 
dipole interactions to be that corresponding to an fee dipole lattice at the given density and 
interpolating the lattice parameter between the LDA values for PbTe and Pb3GeTe4 via 
Vegard's law, we obtain a rough estimate of composition x ~ 0.05, below which there are 
two wells in the potential along the (111) line. Using a different analysis, Katayama and 
Murase estimated the range of validity of their model to be x < 0.1. 

Although our model is classical, it sheds some light on quantum tunneling phenomena in 
the Pbi_^Ge2.Te system at small x by providing a potential energy surface for Ge motion. In 
Figure |Tl|(a), we show the potential energy surface for Ge motion in cubic Pb3GeTe4 at zero 
temperature. There is only a single deep well for the Ge ion and thus no quantum tunneling. 
A more relevant potential energy surface is that for a single Ge impurity in PbTe. While 
our calculations were performed on a structure with 25% Ge, a crude extrapolation to the 
isolated impurity case is possible by (1) setting all Ge-Ge interaction terms equal to zero and 
(2) setting the lattice parameter to the value for pure PbTe. When this is done, the potential 
energy surface for Ge atom motion is that shown in Figure |ll](b). There are now eight wells 
with minima along cubic [111] directions, about 7 meV deep with respect to the centered 
position and with 1 meV barriers between wells. The zero-point vibrational energy of a Ge 
atom in this potential is about 20 meV, so it would not be localized in a single well. At 
some nonzero concentration of Ge impurities, we expect the Ge-Ge off-centering interactions 
to overcome the zero-point vibrational motion leading to an ordered rhombohedral ground 
state. 

The KM model includes not only Ge off-centering, but coupling of Ge motion to the PbTe 
transverse optical phonons. In our calculation of normal modes, the coupling of Ge motion to 
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other ionic motions is given in the dynamical matrix, which is then diagonahzed. Compared 
to the KM model, we both reduce the number of degrees of freedom and incorporate normal 
modes that are more appropriate for a material containing substitution. Consider our Fis 
eigenmode (Table All Pb motion is this eigenmode is in the opposite direction of all 
Te motion, exactly as would be the case of the zone center PbTe TO phonon. The motion 
of one Te atom, however, is much larger in magnitude than all other Pb and Te motions, 
reflecting the importance of localization phenomena. 

It is instructive to compare our results for Pb3GeTe4 with similar results for pure PbTe. 
Linear response calculations on PbTe were carried out using the experimental rocksalt struc- 
ture as the reference structure and the LDA lattice parameter. All normal modes are stable, 
consistent with the experimental fact that PbTe is cubic down to zero temperature.0i It is 
also consistent with the fact that the instabilities in ternary Pb3GeTe4 are associated with 
the Ge ions. 

For pure PbTe, we obtain Born effective charges of +5.84 for Pb and -5.84 for Te and an 
electronic dielectric constant of 33.0. These are in excellent agreement with the experimental 
values of = ±6.0 and eoo = 32. 8@ and suggest that our values for Pb3GeTe4 are also 
reliable. Our calculations show that e^o for Pb3GeTe4 is about 40% higher than for PbTe. 
When Ge is substituted for Pb, the Pb effective charge changes little and remains nearly 
isotropic. The Ge effective charge itself is significantly larger than that for the Pb it replaces. 
The Te(2) Z* also changes substantially and becomes markedly anisotropic. An analogous 
anisotropy has been observed for the oxygen Z* in ABO3 perovskite ferroelectrics. It is 
interesting that the two Z* components of highest magnitude are associated with the two 
most significant ionic displacements in the LWF, namely Ge off-centering along with axial 
motion of the first neighbor Te ions. 

The values of the Born effective charges (and eoo) in Pb3GeTe4 do not remain constant as 
the structure distorts from its reference structure to its ground state. We set the structure 
of Pb3GeTe4 to the ground state of our model {caa = 6.37 x 10^^; Ca/s = 6.03 x 10^^, a 7^ 
P;^ = (0.255,0.255,0.255)) and recalculated the Born effective charge tensors. The results 
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are shown in Table [VI]| . The magnitude of the Born effective charges are generally smaller in 
the distorted state, sometimes much smaller. In particular, Z* for Ge along (111) goes from 
+8.02 to +1.52; the electronic contribution to Z* goes from positive to negative. Previous 
studies have shown similar, but smaller decreases in magnitudes of Z* in BaTiOai^l and 
KNbOgil as the structure is distorted. A common trend is that the atoms with the largest 
relative motion in the structural instability have the largest changes in Z* . The magnitude of 
600 also decreased substantially as the structure distorted. A similar, but relatively smaller, 
decrease has been calculated for KNbOsiil. 

The large difference between for our Pb3GeTe4 model and the experimental Tc for 
the disordered system at the same composition merits discussion. The first source for the 
discrepancy is the approximations used in the first-principles energy calculations. Approxi- 
mations in the ab initio calculations include the use of pseudopotentials, LDA, finite k-point 
sampling, a finite cutoff energy and a finite real-space grid for charge density. The second 
source for the discrepancy is the net effect of the neglected anharmonic effects, which may 
be important in Pb3GeTe4 because the relative displacements of the ions is so large (order 
0.5 A). For example, the dependence of on distortion discussed in the previous para- 
graph will lead to a nonlinear dependence of Pi on ^j, which in turn will lead to anharmonic 
intersite interactions. Such anharmonic effects could either raise or lower and are worth 
further study. The third source for the discrepancy is our use of the LDA lattice parameter 
rather the (unknown) experimental oncH. Typically, ah initio phonon calculations using 
experimental lattice parameters are more accurate than those based on the LDA lattice 
parameters.0 In section 3, we estimated that the LDA lattice parameter a for Pb3GeTe4 is 
about 0.5% smaller than the experimental value, so more accurate results might be obtained 
by raising a. In general, raising the lattice parameter will "soften" the phonons and thus 
raise T^] in the present case, this would increase the discrepancy between our and the 
experimental value. 

Finally, there is the problem of extending the results for an ordered structure to those 
for a disordered alloy. We have, in effect, considered only one term in Eq. |2.1| , which is not 
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in fact "typical". For example, in the cubic configuration no relaxation is possible. Other 
configurations at the same composition do allow relaxation. We have preliminary ab initio 
results that show that the relaxation in disordered Pbi-ajGe^^Te is dominated by a 0.11 to 
0.14 A inward displacement of each Te in the first neighbor shell of each Ge ion, where such 
relaxation does not conflict with the relaxation about a second Ge ion. Such relaxation 
alone lowers the energy per Pb3GeTe4 unit as much as 20 meV per atom with respect to 
the undistorted cubic structure. Locally, the relaxation mimics a decrease of the lattice 
parameter, which should "harden" the unstable Ge-dominated modes and thus decrease Tc. 

The large energy differences between different configurations of Pb3GeTe4 upon relax- 
ation shows the necessity of including different chemical configurations in order to correctly 
model the structural phase transition in the disordered system. This is similar to the case 
of alloy phase diagrams!!, where order- disorder transitions predominate. There, it has been 
shown that lattice entropy effects must be included to obtain correct transition temperatures 
for order-disorder transitions.SJ5i Here, we need to include configurational entropy in order 
to study a structural phase transition. 

The situation for Pbi_a;Gea.Te is complicated by the existence of a miscibility gap, 
with the peak of the exsolution dome at about 840 kI. For X = 0.25, as in this work, 
phase separation to Ge-rich and Pb-rich phases occurs at about 770 although a mixed 
phase of this composition that is metastable at low temperatures can be obtained by rapid 
quenchingS. Thus, the experimental result of ~ 350K for the structural phase transition 
in Pbo.75Geo.25Te actually applies to a metastable compound and the thermodynamic the- 
ory for the disordered compound is complicated by the fact that the equilibrium partition 
function (Eq. |2.1|) does not apply. 

Finally, we discuss the form of the model that should describe disordered Pbi^^Ge^Te. 
For each configuration in Eq. p.l| , it is possible in principle to repeat the procedure used in 
this work to develop a configuration specific model. Based on the Ge off-centering picture, 
we make the following conjectures about the nature of all such models. (1) There is a basis 
for all unstable modes that involves a vector LWF centered on each Ge. (2) The displacement 
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pattern corresponding to each LWF will be strongly localized. (3) The interactions between 
local polar distortions centered on different Ge will favor ferroelectric ordering. 

Given enough models for individual configurations, it should be possible to extract a 
"supermodel" that applies to all Pbi.^jGe^Te configurations. In the same way that energy 
for a disordered alloy can be given in terms of a cluster expansionil, the values of the model 
parameters could also be given in terms of cluster expansions. Work is in progress to develop 
models for other Pbi^^Ge^^Te configurations and to incorporate the effects of relaxation into 
these models. 

VI. CONCLUSIONS 

As a first step towards a first-principles study of the effect of substitutional disorder 
in Pbi_a;Gea;Te, we have developed a model Hamiltonian for an ordered Pb3GeTe4 system. 
We have determined the parameters in this model from ab initio calculations. From a 
classical Monte Carlo simulation, was found to be about 620K. Comparison of the present 
results with those obtained for other configurations will demonstrate how transferable the 
parameters of the model Hamiltonian for this ordered supercell are to other configurations, 
show how much Tc depends on configuration, and allow a model for the disordered system 
to be developed. 

VII. APPENDIX 

The use of a plane wave basis set and periodic boundary conditions in an ab initio 
calculation introduces two kinds of error. The first is due to the finite energy cutoff for the 
plane waves; the second is due to the discretization of wavevectors that is a result of using a 
finite periodic cell in real space. In this paper, we use Pulay corrections to compensate for 
the second source of error. 

The specific method for applying Pulay corrections used here was taken from Rignanese 
et a/.0. In this method, the density functional total energy Etot of a periodic system at 
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fixed volume Vq is measured at several cutoff energies Ecut and treated as a function of the 
average number of plane waves per k point, Npw at those values of Ecut- Rignanese et al. 
suggest using cutoff energies of E^ut (the fixed value for a set of calculations of different cells), 
Ecut — 3% and Ecut + 3% and then fitting through the resulting total energies Etot[Npw, Vo] 
via the following function: 

Etot [Npw, Vol = EZ^ + exp(ao + aiNpw) • (7. 1) 

Finally, the results obtained for the single cell volume Vq are used to determine the Pulay 
corrections to total energy for any cell volume Vi via the correction: 

EUEcut, V,] ~ EfjEcut, V,] + Etot[^Np^{Ecut, Vi),Vo] - Etot^N%w{Ecut. V^),V^\ (7.2) 

Ef^^Ecuti'^-^ is the total energy measured without corrections. A^pvi/ is the corresponding 
average number of plane waves per k point. iVpvF "continuous" number of plane waves 

at Ecut-, determined by multiplying the density of states in reciprocal space by the volume 
of the sphere containing plane waves of energy less than Ecut and is given by 

iVp^(i?cni, V) = ^{2Ecutf'\ (7.3) 

when V and Ecut are in atomic units. From the value of N''p-^{Ecut,Vi) obtained, the 
corrected total energy E'j:^^[Ecut^ ^i] can be computed. 

By measuring the total energy of Pb3GeTe4 at a = 6.375 A and PbTe at a = 
6.275 A at Ecut = 291, 300 and 309 eV, the following fitting forms were obtained: 

PbTe : Etot\Npw, Vq = 259.051^] = -1296.5850 + 6.244 exp(-1.95610-^iVp,^) (7.4) 
Pb3GeTe4 : Etot[Npw, Vq = 247.051^] = -1308.4464 + 8.452 exp(-9.30710~^iVp,^) 

These results were used to correct all the CASTEP 2.1 total energy results used in this 
article. From the expression ( |7.4| ) for Etot, and the value of A^pvi^ = 2911.75 corresponding 
to our calculation for Pb3GeTe4 at a = 6.275 A, we obtained Etot = Uq/N = —1302.1725 eV. 
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FIGURES 

FIG. 1. Crystal structures of the observed phases of pure GeTe. Eight atom regions are shown 

for (a) the high-tcmpcraturc rocksalt structure and (b) the low-temperaturc rhombohcdral struc- 
ture, which is related to (a) by a rhombohedral lattice distortion and a relative displacement of 
the sublattices, shown by arrows. Both distortions are greatly exaggerated for clarity. 

FIG. 2. The high-symmetry reference structure of the ordered cubic Pb3GeTe4 eight-atom unit 
cell, produced by replacing a cubic superlattice of Pb ions (striped circles) in the PbTe rocksalt 
structure with Ge ions (solid circles). 

FIG. 3. The z component of the estimated lattice Wannier function of cubic Pb3GeTe4 (doubled 
for clarity). Only ions with displacements greater than 0.1 times the displacement of the central 
Ge ion are shown. Note that the unit cell outlined is translated with respect to that in Figure 2. 

FIG. 4. (a) Comparison of frozen phonon and fitted results for the unstable F mode subspace 
energy; strain fixed at zero, (b) Contour plot of 110 plane cross section of eighth order expansion 
of unstable F mode subspace energy. The contour interval is 20 meV and the grid spacing is 0.2. 

FIG. 5. Energy vs. strain for (a) uniform dilation [e^x = ^yy = ^zz), (b) tetragonal strain 
i^xx = ^yy = —^zz/^) ^'Hd (c) rhombohedral strain [cxy = exz = ^yz)- 

FIG. 6. Energy change due to local polar distortion of amplitude 0.25 on strained cells. 

(a) uniform dilation; ^ = \^\z, (b) tetragonal strain; ^ = (c) rhombohedral strain; 

FIG. 7. Contour plot of 110 plane cross section of unstable F mode subspace energy, minimized 
with respect to uniform strain. The contour interval and grid spacing are the same as in Figure 4. 

FIG. 8. Result of Monte Carlo simulations, (a) Rhombohedral strain and (b) lattice parameter 
versus temperature. The results are consistent with a second order phase transition with near 
620 K. 
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FIG. 9. Phonon dispersion relations for the cubic Pb3GeTe4 model Hamiltonian. 



FIG. 10. Distribution of ^ in Monte Carlo simulation, (a) Contour plot in 110 plane at Tc+SOK. 
(b) Same as (a), at T = T^-SOK. (c) First neighbor Ge-Te distribution function at T = Tc+SOK. 
(d) Same as (c) for T = Tc-SOK. Compare (a) with figure 4(b). The eight [111] directions are 
preferred for local distortions, even above Tc, in agreement with the Ge off-centering picture. 

FIG. 11. (a) Potential energy surface seen by single Ge ion at T = in Pb3GeTe4 model 
(contour interval 50 meV). (b) Potential energy surface seen by isolated Ge impurity in PbTe, 
based on extrapolation of our model (contour interval 5 meV). 
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TABLES 

TABLE I. Ab initio frequencies (in cm~^) of selected normal modes in cubic Pb3GeTe4. The 

frequencies for optical Fis modes are for the transverse modes. 



Symmetry Label: 


ri5 


r25 


Xy 


^5' 


Mil 


M5/ 


Rl5 




648 i 


176 


478 


593 i 


547 i 


395 


429 







630 


518 


163 


315 


603 


437 




195 




721 


308 




688 


683 




292 






622 




736 






518 
















662 















TABLE II. Displacements (in A) associated with all symmetry-independent modes incorpo- 
rated into the model Hamiltonian. The periodicity of each mode listed except Lis is longer than 
the primitive cell and involves the opposite motion of some ions in neighboring cells. 



Mode 


Pis 




X5/ 


Mr 


M5' 


Rl5 


uj (cm~^) 


648 i 


478 


593 i 


547 i 


395 


429 


q (27r/a) 


(0,0,0) 


(0,0,0.5) 


(0,0.5,0) 


(0.5,0.5,0) 


(0,0.5,0.5) 


(0.5,0.5,0.5) 


u (Ge(0,0,0)) 


0.8878 z 


0.8258 z 


0.8835 z 


0.8654 z 


0.8207 z 


0.8578 z 


u (Te(0,0,0.5)) 


-0.4467Z 





-0.4668i 


-O.SOlOi 








u (Te(0,0.5,0)) 


-0.0575i 


0.0680 z 














u (Te(0.5,0,0)) 


-0.0575Z 


0.0680 z 


-0.0318Z 





0.0555Z 





u (Pb(0.5,0.5,0)) 


0.0156 z 


-0.5556i 














u (Pb(0.5,0,0.5)) 


0.0284 z 





0.0237 z 








-0.3569X 


u (Pb(0,0.5,0.5) 


0.0284 z 











0.5684y 


-0.3569y 


u (Te(0.5,0.5,0.5)) 


-0.0611Z 











-0.0163y 
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TABLE III. Ionic displacement pattern corresponding to z component of estimated lattice 
Wannier function Wzji for c-Pb3GeTe4. The lattice Wannier function transforms according to the 

Fis representation for point group Oh with center R. Orbits of equivalent atoms for Wzn. ai's labeled 
according to the International Tables^^ convention for point group A/mmm. Components whose 



valuers 


are vanish by symmetry. 


while tliosc^ whose values aix 


^ 0.0 are cnily zero by ai")|")r(Lxiiuati(ni. 


Atom 


r j — R (units of a) 


Multiplicity (orbit) 


u(A) 


Ge 


(0,0,0) 


1 


(0,0,0.8557) 


Te 


(0,0,0.5) 


2(a) 


(0,0,-0.2352) 


Ge 


(0,0,1) 


2(a) 


(0,0,0.0122) 


Te 


(0.5,0,0) 


4(c) 


(0,0,0.0043) 


Pb 


(0.5,0,0.5) 


8(f) 


(0.1165,0,0.0065) 


Te 


(0.5,0,1) 


8(f) 


(0.0,0,-0.0133) 


Pb 


(0.5,0.5,0) 


4(b) 


(0, 0, -0.0675) 


Te 


(0.5,0.5,0.5) 


8(e) 


(-0.0020, -0.0020, -0.0076) 


Pb 


(0.5,0.5,1) 


8(e) 


(0.0,0.0,0.0357) 


Ge 


(1,0,0) 


4(c) 


(0, 0, -0.0006) 


Te 


(1,0,0.5) 


8(f) 


(0.0,0,0.0034) 


Ge 


(1,0,1) 


8(f) 


(0.0,0,0.0017) 


Te 


(1,0.5,0) 


8(d) 


(0, 0, -0.0008) 


Pb 


(1,0.5,0.5) 


16(g) 


(0.0,0.0128,0.0003) 


Te 


(1,0.5,1) 


16(g) 


(0.0,0.0 - 0.0012) 


Ge 


(1,1,0) 


4(b) 


(0,0,0.0009) 


Te 


(1,1,0.5) 


8(e) 


(0.0,0.0,-0.0004) 


Ge 


(1,1,1) 


8(e) 


(0.0,0.0,-0.0009) 
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TABLE IV. Born effective cliarges for cubic Pb3GeTe4. For Ge and Te(l), ^* is isotropic, 
while for Pb and Te(2), the parallel and perpendicular components refer to the local tetragonal 
axis. 



Species 


Z* 


zi 


Pb 


+6.15 


+6.06 


Ge 


+8.02 


+8.02 


Te(l) 


-5.82 


-5.82 


Te(2) 


-7.73 


-6.37 



TABLE V. Harmonic energies of effective Hamiltonian modes (|^| = 1 on each Ge site). Units 
are eV per unit cell. Displacements of ions are given in Table II. The energy for Fis was determined 
via a frozen phonon calculation, the others are determined via nonzero q linear response. 



Mode 






\J sr 


Energy via Eqs. 3.7-3.9 


ri5 


-1.652 


-0.383 


-1.269 


A + bL + 2 bT + 2 (di + dTi) + 2 dT2 + 4 (fL + 2 fr)/3 


Xr 


1.220 


0.885 


0.335 


A - bi + 2 br - 2 (di + dTi) + 2 dr2 - 4 (fi + 2 fT)/3 


X^' 


-1.387 


-0.443 


-0.944 


A + bL - 2 dT2 - 4 (fi + 2 fT)/3 


Mi> 


-1.206 


-0.489 


-0.717 


A - bi - 2 dT2 + 4 (fi + 2 fT)/3 


Mr, 


0.811 


0.245 


0.599 


A + 1)L 2 hr 2 (fU + d:,-i) + 2 dr2 + 1 (.fi + 2 rr)/3 


Rib 


0.927 


0.000 


0.927 


A - bi - 2 bT + 2 (di + dTi) + 2 dT2 - 4 (fL + 2 fr)/3 
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TABLE VI. Parameters in the effective Hamiltonian for Pb3GeTe4 (units of eV per unit cell). 



A 


-0.1765 


hL 


-0.7922 




215.8 


B 


4.501 


bT 


-0.2148 


Cl2 


9.374 


C 


6.165 




0.0766 


C44 


108.1 


u 


— / .ZOo 






50 


— o.oU 


E 


4.239 


dT2 


-2.98 X 10-2 


91 


-12.4 


Z* 


12.10 


h + 2 It 


2.55 X 10-2 


92 


-9.99 




46.7 


a 


6.275 A 







TABLE VII. Born effective charges and eoo for ground state of cubic Pb3GeTe4. Atomic 
positions are given to the nearest tenth. The principal directions are given by the subscripts. They 
are exact by symmetry for Ge, Te(l) and eoo, but only approximate for Pb and for Te(2). The 
corresponding values for the high symmetry reference structure are given in Table IV and shown 
here in parentheses. 



Pb (0.0,0.0,0.5) Z* 


+5.86001 (6.15) 


+5.84ijo (6.06) 


+4.98110 (6.06) 


Ge (0.0,0.0,0.0) Z* 


+1.52iii (8.02) 


+3.89ijo (8.02) 


+3.89ii2 (8-02) 


Te(l) (0.5,0.5,0.5) Z* 


-5.67iii (-5.82) 


-5.82^^0 (-5.82) 


-5.82^12 (-5.82) 


Te(2) (0.5,0.5,0.0) Z* 


-3.65ooi (-7.73) 


-5.OO1T0 (-6-37) 


-5.42iio (-6.37) 


^00 


3I.O111 (46.7) 


30.9ijo (46.7) 


30.9i,2 (46.7) 
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Pb 

Ge 
OTe 





-600 i = 



r X M R 





(b) 



